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The process of gap formation by a growing planetary embryo embedded in a planetes- 
imal disk is considered. It is shown that there exists a single parameter characterizing 
this process, which represents the competition between the gravitational influence of the 
embryo and planetesimal-planetesimal scattering. For realistic assumptions about the 
CN| ■ properties of the planetesimal disk and the planetary embryo, a gap is opened long be- 

fore the embryo can accrete all the bodies within its region of influence. The implication 
■r^- | of this result is that the embryo stops growing and, thus, large bodies formed during 

^ \ the coagulation stage should be less massive than is usually assumed. For conditions 

expected at 1 AU in the solar protoplanetary disk, gap formation is expected to occur 
around bodies of mass < 10 24 g. The effect of protoplanetary radial migration is also 
discussed. 
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1. Introduction. 

The formation of planets is one of the most complex problems in astrophysics, involving accu- 
mulation of bodies over some 45 orders of magnitude in mass — from dust grains to giant planets. 

One issue which has received a lot of attention is the formation of planetary embryos by the 
accretion of planetesimals. From the perspective of dynamics we call an object an embryo when 
it becomes so massive that one can no longer describe its behavior by means of simple kinetic 
theory (in this paper we use names embryo, protoplanet and massive body interchangeably). In 
other words, in the presence of embryos the multiparticle distribution function cannot be taken as 
a product of one-particle distribution functions; the gravitational influence of the embryo is strong 
enough to affect the distribution of planetesimals with which it interacts. For example, a gap could 
form around the embryo. To study properties of systems containing embryos one must either resort 
to N-body simulations or try to account properly for their influence on the underlying planetesimal 
population and on each other. 
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In the standard scenario, protoplanets grow in orderly (Safronov, 1972) or runaway fashion 
(Wetherill & Stewart, 1989; Wetherill & Stewart, 1993) by accreting planetesimals from the pro- 
toplanetary nebula. After the largest bodies become embryos, they open gaps and accretion slows 
or stops. Thereafter these massive bodies evolve more slowly as gravitational encounters perturb 
them into crossing orbits and violent impacts occur, thus gradually forming more massive objects. 
It is a common belief now that the Earth-type planets and rocky cores of the giant planets were 
formed by this two-stage process. 

The question which is not very often addressed is where the boundary between these two stages 
occurs. This is an important issue, because the answer tells us the final mass of the objects which 
further evolve through chaotic collisional evolution, as well as the number of such bodies, and the 
conditions for which statistical treatments of collisional evolution are valid. 

The standard paradigm for determining the embryo mass (Lissauer, 1987; Weidenschilling et 
al. 1997) presumes that planetary growth stops when the embryo "eats up" all the planetesimals 
within a "feeding zone" - an annulus of radial width of one Hill radius. Here the Hill radius is 
defined as 

with a being the distance between the massive body with mass M p and the central star with mass 
M c . If the surface mass density of planetesimals is So, then this "isolation mass" is M p = M{ s ~ 
27rai?f/£o, which means that 

MtS ^F^ = 1Xl ° g 1.2 x 10- M J UeJ • (2) 

The estimate of Y>qo? is made for 1 AU and is based on standard assumptions for the protosolar 
nebula: So ~ 20(a/l AU) -3 / 2 g cm -2 , implying that ~ 1 per cent of the minimum mass Solar 
nebula is contained in solids (Hayashi 1981). 

The usual assumption is that before the mass of the largest protoplanet reaches Mj S , the 
distribution of planetesimals is basically homogeneous. In some cases this is not a reasonable 
assumption. In particular Ida & Makino (1993) demonstrated using N-body simulations that a 
protoplanet could scatter planetesimals strongly if it is massive enough and, thus, clear a zone 
around it which is free from any solid bodies. The mass required to clear a gap in this way is not 
simply related to Mj s . Kokubo & Ida (1998) later emphasized the importance of rapid heating 
of the planetesimal population by the forming planet in slowing down the subsequent accretion of 
planetesimals. 

The process of clearing a gap in the planetesimal disk around a massive body is analogous to 
gap formation in gaseous disks (Takeuchi et al. 1996), which results from a competition between 
viscous spreading of the disk and gravitational interactions with the protoplanet. In a planetesimal 
disk the role of viscosity is played by mutual scattering of planetesimals. One can easily estimate the 
planetary embryo mass determined by this process. Let us assume that the planetesimal random 
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velocities are small - the disk is cold. Let Q = (GM c /a?) l l 2 be the orbital angular velocity. If tuq is 
the mass of each of two planetesimal and rjj = a(2mo/M c ) 1 / 3 is the corresponding Hill radius, then 
the typical displacement in a close encounter of these planetesimals on circular orbits separated by 
< th is ~ rjj, and the typical random velocity kick is f2r#. Similarly, planetesimals within Rh 
from the massive object get kicked by ~ Rh with frequency ~ Q(Rh/cl) (the inverse of the synodic 
period). If they are not able to diffuse back this distance during the time interval between kicks by 
the embryo, then a gap forms. 

The "viscous" spreading distance is ~ rnK 1 / 2 during one synodic period, where K is the 
number of collisions of a given small body with other planetesimals between consecutive approaches 
to the massive body. Assuming that the thickness of the planetesimal disk is ~ ru , one can easily see 
that K ~ r H (Y, /m )(a/R H ) = a 2 S /(m M p M c ) 1 / 3 . This implies that a gap in the planetesimal 
disk opens when R H > r 2 H K, or when 



M P M C 1/3 



> 1, (3) 



f(v/Q,r H )^oa 2 ml /3 

or, alternatively, when 

/ \ 1/3 

M p > M crit ~ f(v/nr H )Z a 2 . (4) 

Here /(f/f2r#) is a dimensionless function characterizing the effect of the planetesimal velocity 
dispersion v; we expect f(x) ~ 1 for x < 1 and f{x) ~ x~ 2 \nx for i> 1 (see §4.2). If the mass 
given by equation (4) is smaller than the isolation mass given by equation (2), the accretion stops 
because the protoplanet forms a gap, rather than because it consumes all the bodies in its feeding 
zone. Assuming typical values for protoplanetary disks one can get 

It is obvious from this estimate that gap formation could be very important in slowing down 
planetary accretion. 

In this paper we analytically study the process of clearing a gap around a massive body in 
a planetesimal disk. We use an approach to treating the surface density evolution that was first 
developed by Petit & Henon in their seminal series of papers (1987a, 1987b, hereafter PH, 1988). 
In §2 and Appendix A we derive a generalized form of their evolution equation, including the fluxes 
produced by the protoplanet and those generated by mutual gravitational perturbations between 
the planetesimals of the swarm. 

In §3 we describe the solutions of the evolution equations for cold planetesimal disks, and com- 
pare our results with those obtained using N-body simulations. We comment on the applicability 
of our findings to the planet formation in the early Solar System, and describe briefly the relation 
between the surface density and planetesimal velocity dispersion evolution in §4. 
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2. Derivation of the general equation. 



All the following calculations assume a Keplerian disk, although they could be easily extended 
to the case of an arbitrary rotation law. 

We consider a disk of bodies (we will refer to them as planetesimals, but these could be other 
bodies, such as planetary ring particles) with N(m,r,t)dm = Y±(m,r,t)dm/m being the surface 
number density of particles with mass between m and m + dm, whose guiding centers move at 
a distance r from the central body. It is important to keep in mind that r is the guiding center 
radius rather than the instantaneous radius. The instantaneous surface number density can only 
be obtained if T,(m,r,t) is supplemented by the random velocity distribution of planetesimals. 

We also assume that a single massive body with mass M p moves on a circular orbit in this 
planetesimal swarm (we take orbit to have a fixed radius and we comment later on the effects of 
migration) and we assume that its mass is much larger than the masses of the individual swarm 
particles. The mass of the central body is M c and the distance of the planet from the central body 
is a. We will also use the relative masses of the bodies with respect to M c : \i v = M p /M c for the 
planet and \i = m/M c for the planetesimals with mass m. 

The interactions between particles in the gravitational field of a central body are described by 
Hill's equations, which are valid in the limit fi,n p <C 1, which is always true in problems which we 
will study. It was demonstrated by Henon & Petit (1986) that in this case the motion of nearby 
gravitationally interacting particles can be separated into center-of-mass motion, which is invariant 
during the interaction, and relative motion. If one defines new dimensionless coordinates where 
all the distances and relative velocities are normalized by a(pL\ + fi^ 1 ^ 3 , then the equations of 
relative motion of particles 1 and 2 do not depend on their masses in these coordinates. Let us 
set h to be the distance between the guiding centers of interacting particles in these coordinates 
and P(h, Ah)dAh to be the probability of having a change in h in the range (Ah, Ah + dAh in 
an encounter. Then P(h, Ah) does not depend upon the masses of particles involved in a collision, 
but does depend on the random velocity distribution function of the planetesimals. 

In Appendix A we derive the general equation of the surface density evolution, which in 
many aspects parallels the derivation of equation (44) in PH. Let us set Ni = N(m,i,r,t). Let 
A = (r/2)(dfl/dr) be the function determining the local shear, A = —(3/4)0 for a Keplerian 
rotation law. Then the surface density evolution is given by 



oo 




oo 



azvi(r) 

dt 



2\A\a 2 / dm 2 ( f i 1 + fi 2 f/ 3 / dh\h\{ N^N^r - ( W + fi 2 ) 1/3 rh] 



o 



— oo 



oo 



d{Ah)P{h, Ah)Ni [r + D(Ah)]N 2 [r + D(Ah) - (m + fi 2 ) 1/3 ah] 



(6) 



— oo 
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where 

D(Ah) = - f° 3 Aft. (7) 
(Ail + ^2j 2/d 

Note that the factor r is replaced with a in equations (6), (7), where it is appropriate, because at 
this level of approximation there is no difference between them, since they are both much larger 
than the Hill radius. 

This differs from the equation derived in PH because it does not assume that surface density 
varies slowly on scales of the order of the Hill radius. If we made this assumption and expanded 
N(m, r, t) up to the second order in h locally we would reduce equation (6) to the one derived in 
PH, equation (44), which we reproduce here: 



oo 

= Wf** [ 2W ,, W 3 | (n^) + / 2 ^f 5575 



dNi 





li\ d 2 (N 1 N 2 ) 
dr 2 



(8) 



Here I\ and I2 are dimensionless moments of the probability distribution P(h,Ah): 

00 00 00 

h = j \h\hdh j d(Ah)AhP(h,Ah) = J \h\hdh(Ah), (9) 

—00 —00 —00 

00 00 00 

h= J \h\dh j d(Ah)(Ah) 2 P(h,Ah) = J \h\dh((Ah) 2 ), (10) 

—00 —00 —00 

and symmetry dictates that 

00 00 

J \h\dh J d(Ah)AhP(h, Ah) = 0. (11) 

—00 —00 

For a cold disk it was demonstrated by PH that 

l x = -3.07, I 2 = 17.72. (12) 

We can now easily include the effect of a massive body on the surface density evolution. To 
do this we take surface density to consist of two parts: one representing a continuous distribution 
of small masses, corresponding to planetesimals, and another arising from the massive body. One 
can write the contribution from the embryo in the following form: 

N em (m,r,t) = ^-5(m-M p )5(r-a). (13) 
We neglect migration of the embryo, so we assume a is fixed. 
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Substituting (13) into (6) and assuming that the embryo is much more massive than any of 
the planetesimals, M p S> m, we get 



-2Aa 2 J dm 2 (tn+to) 2/3 J dh\h\ \ N 1 (r)N 2 [r - (m + fi 2 ) 1/3 ah] 

-oo ^ 

oo 

- J d(A/t)P(/i,A/i)iV 1 [r + D(A/i)]^ 2 [r + L>(A/i) - ( m + fi 2 ) 1/3 ah] } (J I.) 



A 



na 



Ni(r)\r-a\ - 



OO / 

-j^ / dnJVi(ri)|ri -a\P\ — — f 



H 

/3 



We can make other simplifications taking the following into account. In our particular problem 
the relevant length scale for any structure is the Hill radius of the massive body Rh- Planetesimals 
can get kicks when interacting with the massive body which change their guiding centers by ~ Rh- 
At the same time mutual interactions between the planetesimals are unable to produce such large 
displacements (since /x p S> /ij). One can thus assume that for planetesimal-planetesimal interactions 
the surface density varies only slowly on the scale ru <C Rh, and locally expand the first part of 
the r.h.s. of equation (6) in a Taylor series in h, as was done in equation (8). At the same time the 
second part of the r.h.s., representing the interaction with the large body, cannot be simplified in 
a similar way. 

We also make some additional changes: we move the origin of r to a (simply set r — a = r 1 ) 
and switch from r' to a dimensionless distance from the planetary embryo H = r'/(^ 3 a). Then 
we get 



/4 /3 dNi(H) 
Aa 2 dt 



/ 



dm 2 



, l7 o d ( dN 2 \ ii\ d 2 (N 1 N 2 ) 



7ra 2 



oo 

N!(H)\H\- J dm Ni(Hi)\Hi\P(Hi, H — H±) 



(15) 



In deriving this form of the evolution equation we only assumed that /i p S> So, this 

nonlinear integro-differential equation can adequately describe the evolution of the surface density 
of planetesimals in the disk-protoplanet system. 



2.1. Single mass planetesimals. 



The constituent bodies of planetesimal disks are likely to have quite a broad range of masses. 
However, right now we are going to concentrate on the simple case of single mass planetesimals, that 
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is we assume all planetesimals to have a unique mass mo. Then N(m, H) = a(H)(T,Q/mo)5(m—mo), 
where So is the surface mass density of particles at infinity which we take to be a reference value 
(it follows then that cx(oo) = 1). Substituting this assumption into equation (15) and performing 
an integral over one obtains that 



1 da 8 2 a 2 

A 



oo 

a{H)\H\- J dH\ a(Hi)\Hi\P(Hi, H — Hi) 



(16) 



I dr dH 2 
where /j,q = mo/M c and 

/ = /! + — = - f \h\(2hAh + (Ah) 2 )dh. (17) 
Z z J 

— oo 

The new time variable r is defined as 



oo 



r = l j W here t = , (18) 



and [cf. eq. (3)] 



M p 



(19) 



2 1 /3 vrSo a 2 /A/o /3 ' 

One should notice that the first term in the r.h.s. of equation (16) and the expression in brackets 
are both dimensionless; all of the dimensional information is hidden in A and r. 

For a cold disk (rms velocity dispersion of planetesimals in r-direction v r fJr#) we obtain 
using (12) that 

/ = 5.79, (20) 

and, thus, 

A = 0.0436 — ^ytj, v r < ttr H . (21) 
For a hot disk {v r S> f2r#) one has [see the discussion after equation (33) in §4.2] 

Q2 r 2 / 2 \ 3/2 



7 = 29.8 ^ In A, with A ~ j . (22) 

It is assumed here that the ratio of vertical to radial random velocity dispersions in a planetesimal 
disk is equal to 0.5. Using (22) one finds that 

M v 2 

A = °-° 085 v 2 P i/3 2 ;2 ( ln ^ , v r » nr H . (23) 

In the intermediate regime, when v r ~ Qrn, there is no analytic expression for I and one has to 
interpolate between the two asymptotic behaviors given by (20) and (22). 
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The parameter A quantifies the influence of the planetary perturbations on the uniformity of the 
planetesimal disk. The first term on the r.h.s. of (16) describes the nonlinear diffusion of particles 
due to mutual gravitational scattering, and tends to iron out any initial inhomogeneities. The 
second term represents the effect of the planet, which tends to carve out a gap in the distribution 
of planetesimals. The steady state originates when these two effects balance each other. 

When the protoplanetary mass is large, A is also large and the expression in brackets dominates 
the evolution. It leads to gap formation. On the contrary, if the planetary mass is small we can 
neglect the corresponding term in the r.h.s. of (16) and obtain a nonlinear diffusion equation, which 
drives the planetesimal distribution towards a homogeneous state. So, one can say that a gap (or 
at least a significant depression in the surface density of the planetesimals) is formed when A > 1. 
From equations (16) - (19) one can derive the characteristic time required for a gap to form when 
A>1: 

where T = 2ir/£l is the orbital period of the embryo. Note that t open is approximately the synodic 
period of a body at Rh from the embryo and, thus, is independent of the planetesimal mass mo, 
the surface density S, and the numerical factor /. 

In §3 we confirm these arguments by solving equation (16) numerically. 



3. Numerical results. 

3.1. Solution of the equation of evolution. 

We solved equation (16) in a simplified setting, in which we completely neglect the velocity 
evolution of the planetesimal population. Thus, we assume the integral / embodying all the kinetic 
properties of planetesimals to be fixed in time. For simplicity we set / = 1; this only affects the 
timescale of the gap formation by a constant factor and does not change the evolution at all. 

We also assume that planetesimals have very small random motion on the scale of the embryo's 
Hill radius: v <C flRn- This simplifies our treatment a lot, because in this case scattering is 
deterministic so that P(h,Ah) = 5(Ah — h'(h) + h). The behavior of the function h'(h) which 
gives the final semimajor axis difference as a function of the initial difference was described in 
detail by PH. For the sake of convenience we reproduce this dependence in Appendix B. Also, in 
this approximation, the instantaneous surface number density is given simply by cr(H), because 
the guiding center and instantaneous radii coincide. The assumption of a cold disk might be 
reasonable in cases such as the early stages of gap clearing in a planetesimal disk, when scattering 
by the embryo could be in the shear-dominated regime, or all the time in dense planetary rings 
(Petit & Henon 1987a). 

We solved equation (16) with periodic boundary conditions, simply assuming that da/dH = 
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Fig. 1. — The time evolution of the surface density in a cold planetesimal disk with a single massive 
body at H = 0, for two values of the parameter A defined in equation (19). In the top panel the 
case A = 1 is described, in the bottom panel we consider A = 100. The dimensionless time r is 
indicated on the panels; larger r corresponds to a deeper gap. Notice the presence of a bump at 
the center of a forming gap which is due to particles in horseshoe orbits near the massive body. 
The increase in asymptotic surface density at late r is an artifact of the use of periodic boundary 
conditions at H = ±20, which forces J^o °~(H)dH to be conserved. 
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at 77 = ±7. The functional form of P(h,Ah) is given by equation (Bl). We usually assume 
L = 20.0 here (in units of the Hill radius of planet) and take the initial surface number density to 
be constant: a(H,0) = 1. 

In Figure 1 we show the evolution of the surface density with time [we use the dimensionless 
time r given by equation (18)] for A = 1 and 100. In the first case a gap is never actually formed 
and in the steady state there is only a density depression around the planet. Thus the particles 
can be still accreted by the protoplanet, but the efficiency of this process is reduced. 

In the case A = 100 the gap is formed very quickly, which is in general agreement with the 
expected timescale for gap formation (r open ~ A -1 in this case), although it takes some time after 
that for the density distribution to settle to a steady state. 

In both cases one should notice a bump inside the gap which decays with time. It corresponds 
to the horseshoe orbits in the immediate vicinity of an embryo. This is in agreement with Monte- 
Carlo and N-body simulations performed earlier (Petit & Henon, 1988; Tanaka & Ida, 1997; Spahn 
k Sremcevic 2000). 

In Figure 2 we show the final state of the surface density for several values of A, so that one 
can see that gap gets deeper and wider as A increases. It also looks like the condition A = 1 is a 
reasonable approximate criterion for gap formation. 

3.2. Comparison with N-body simulations. 

We may compare our analytical results from equation (19) with N-body simulations by Ida & 
Makino (1993, hereafter IM93) and Tanaka & Ida (1997, hereafter TI97). 

Figure 3 of IM93 and Figure 1 of TI97 show well developed gaps in a gas-free planetesimal 
disk. The gaps seen in N-body simulations are never clean because random motion of planetesimals 
is naturally included, and this permits some of them to be present in the gap. The parameters 
used in the production of these Figures correspond to A « 25// in the first case while in the second 
A « 60/7. In both cases the velocity dispersion of planetesimals is large (v/Urn ~ 20 — 70) and we 
expect 7 <C 1 (see §4.2), meaning that A S> 1 and the condition for gap formation in the distribution 
of guiding centers of planetesimals should be fulfilled. It was also demonstrated by TI that gas 
drag could clear the gap of residual high-velocity planetesimals and thus stop accretion completely. 

In Figure 6 of IM93 there is shown a sequence of scenarios for different ratios of the M p /mo 
- planet to planetesimal masses. In the case M p /mo = 10, when the gap is barely seen at all, 
A i=3 2.2/7; since planetesimals are not strongly heated, presumably A < 1. For M p /mo = 30, 
when the gap becomes pronounced, A 6.6/7 and A ~ 1. Finally for M p /rriQ = 100, when the 
gap is quite significant, A P3 22/7; heating also becomes important and A> 1. In addition, these 
results confirm our prediction below in §4.2 that for A ~ 1, when gap starts to form, planetary 
perturbations begin to dominate planetesimal random velocity stirring within ~ Rh of the planet. 
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Fig. 2. — The final distribution of the surface density for a cold disk with a single massive body, 
for several values of the parameter A: 0.2; 0.5; 1; 10; 100. The higher values of A correspond to 
progressively deeper gaps in the Figure. The increase in asymptotic surface density at late r is an 
artifact of the use of periodic boundary conditions. 
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This comparison shows that estimates of gap formation based on the parameter A are qualita- 
tively consistent with N-body simulations. 

4. Discussion 
4.1. Applications. 

We have established the usefulness of the parameter A in describing gap formation. Now we 
are going to use it to determine the mass of the body which could open a gap in a planetesimal 
disk. This is simply done by setting A = 1. We rewrite this condition as 

M crit = 2 1 /3^ /Soa 2 ; (25) 

where M cr a is the planet mass for which A = 1. One can easily see that this value of critical 
planetary mass coincides with our simple estimate in equation (4) if we set f(v/Qrn) = 2 1 / 3 ttI. 

As we mentioned in §1, it is usually assumed that accretion stops when planet hoovers up a 
zone in the planetesimal disk with width equal to its Hill radius [eq. (2)]. Then the ratio of the 
two critical masses is 

Merit _J_ ( m A 1/3 / M c A 1/6 

M is 27/6^1/2 \^ oa 2) \^ a 2) ■ \ ^ 

We expect the total mass of the circumstellar disk to be of order 0.1 — 0.01 of the mass of the 
central star M c (Osterloh & Beckwith 1995; Mannings & Sargent 2000). For the protosolar nebula 
it is often assumed that surface density of gas at 1 AU is ~ 2000 g cm -2 (Hayashi 1981). The 
mass fraction of heavy elements, which contribute to solid body formation in this disk, is ~ 0.01. 
Using this information we can estimate that M c /(£o« 2 ) ~ 10 5 — 10 6 at 1 AU and the corresponding 
factor does not contribute a lot to the ratio in (26) (it varies roughly from 7 to 10). The factor 
2~ 7 / 6 7r~ 1 / 2 / ss 1.5 for a cold disk [see eq. (20)], and is significantly smaller for hot disks [it is likely 
to be ~ (m /Mp) 2 / 3 , see §4.2]. 

If we take planetesimals to be rocky bodies with radius ~ 50 km and mass ~ 10 21 g, then 
m /(E a 2 ) ~ 10~ 6 - 10" 7 . In the end we obtain that 

M crit /M ls ~ (10- 1 - HT 2 ). (27) 

This result means that accretion is slowed down long before the clearing of the feeding zone. We 
conclude that the masses of planetary embryos are ~ 10 to 100 times smaller than predicted by 
arguments based on clearing the feeding zone [see eq. (2)]. 
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4.2. Random motions of planetesimals 

Equation (16) fully describes the surface density evolution (or its steady-state structure) only 
if it is supplied with information about random motions of particles within the planetesimal disk, 
which determine P(h, Ah) and its moment /. 

For a cold disk (zero velocity dispersion) I is given by (20). In our case the planetesimal 
swarm is unlikely to be cold, because planetesimals will be scattered by the planet, and also will 
scatter each other. An important point to note here is that these two types of scattering probably 
operate in quite different regimes. Indeed, the natural parameter determining the heating regime 
is the ratio of the velocity dispersion to the shear across the Hill radius. For the scattering by 
the planet this parameter is S ~ v/QRh, with Rh = a(M p /M c ) 1 / 3 while for the interaction with 
other planetesimals it is s ~ v/ttrn, with th = o(2mo/M c ) 1 / 3 , so that s ~ S(M p /mo) 1 ^ 3 ■ Thus, 
scattering by the planetary embryo could be in a shear-dominated regime, while mutual scattering 
of planetesimals is practically always in a dispersion-dominated one. 

During a passage within a Hill radius from a massive body, a planetesimal initially on a circular 
orbit gets a significant kick, so that its velocity dispersion increases by ~ QRh- Thus, the part of 
the disk within a Hill radius from the planetary embryo will be heated to v ~ QRh corresponding 
to S ~ 1 in time ~ Q~ 1 (cl/Rh) = f2 _1 /Xp 1//3 . Thus the stirring rate by planetary scattering is 

/ M 3 \ V3 

^(^fe ' fOT S ^ < M > 



dt 



pi \' o 



within a radial distance Rh from the embryo. The same kind of estimate could be derived for the 
excitation of the vertical random motions. 

This heating rate quickly leads to s 3> 1 and, thus, the self-heating of the planetesimal popula- 
tion should be calculated in a dispersion-dominated regime. In particular the integral / in equation 
(16) describing the scattering of planetesimals by their mutual interactions should be calculated in 
this approximation. 

Stewart & Ida (2000) considered velocity stirring in the dispersion-dominated regime. Their 
results demonstrate that the horizontal stirring rate is 

dt 



= n^r H (P vs ), (29) 

self m 



if the vertical velocity dispersion is of the same order as a horizontal one (as we normally expect), 
and they provide a closed analytic form for the stirring coefficient {Pvs) as a function of planetesimal 
velocity dispersion. The coefficient (Pvs) is defined to represent the change in the square of the 
relative eccentricity [see Henon & Petit (1986)], averaged over the vertical and horizontal orbital 
phases, and the velocity distribution of planetesimals (Ida, 1990; Stewart & Ida, 2000): 

/' 3 dTcluj 

Ael(e u ,i u , h, r, u)f(e u , i u )de 2 u di 2 u - \h\dh^^- , (30) 
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where e u ,i u ,T,tu are the relative eccentricity, inclination, horizontal and vertical orbital phases, h 
is the separation of semimajor axes of the interacting particles and f(e u ,i u ) is the corresponding 
distribution function. We assume that e u , i u , and h are properly normalized by the Hill radius for 
particles participating in the collision [in this case our definition of (Pvs) differs from Stewart & 
Ida (2000) by a factor of (m 1 + m 2 ) 4 / 3 /(3M c ) 4 / 3 , where mi and 777-2 are the masses of interacting 
planetesimals] . 

Vertical stirring is described by a formula similar to (29) but with a different stirring coefficient 

/3 cItcIuj 
A^(e u , i u , h, r, iv)f(e u ,i u )deldil-\h\dh^^-. (31) 

Now we can say more about the relation of the integral / to the kinetic properties of the 
planetesimal population. One can easily see that the averaging over all possible Ah for a given 
initial h used in definition (17) is equivalent to averaging over the vertical and horizontal orbital 
phases at a given relative eccentricity and inclination and then over the distribution of relative 
eccentricities and inclinations, that is 

00 

I = \ J Wh J f(e u ,z u )deldi 2 u ^[2hAh + (Ah) 2 }. (32) 

—00 

From the conservation of Jacobi constant one has 2hAh + (A/7) 2 = (4/3)[A(e 2 ) + A(i 2 t )]. Then, 
substituting into (32) we get that 

I=l((Pvs) + (Qvs)). (33) 

This is an important result, because it relates the evolution of the surface density of a spatially 
inhomogeneous planetesimal population to the viscous stirring in a homogeneous disk. 

Using these expressions we can determine the relative role of self-heating of the disk and 
planetary heating in the vicinity of the massive body. Of course, the former dominates beyond 
several Rh from planet, but within ~ Rh of the embryo one can get from equations (19), (28), 
(29), and (33) the simple result that 

{ds 2 /dt)\ pl x 



(ds 2 /dt)\ self 

This means that for A > 1, when a gap starts to form, embryo perturbations begin to dominate 
planetesimal heating within ~ Rh of the planet, i.e. when the planet dominates the surface density 
evolution it also dominates the heating. 

The evolution of kinetic properties of planetesimals could be neglected if there is an effective 
velocity damping due to inelastic collisions between bodies, as in the case of planetary rings (Petit 
& Henon, 1987a), or if there is a strong gas drag. In the planetesimal case, however, gravitational 
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stirring dominates over damping (Kenyon & Luu 1998) and then the fact that a gap in the distri- 
bution of guiding centers is opened does not automatically mean that planetesimals cannot reach 
the planet, because in the course of scattering their velocity dispersion grows as well [see equa- 
tion (34)]. However, the accretion rate will drop anyway at least because of the less pronounced 
focussing (Safronov, 1972; Dones & Tremaine, 1996). Also, inhomogeneous random velocity evo- 
lution or gas drag could remove the residual planetesimals from the forming gap, thus bringing 
their surface density around the embryo to zero and shutting down accretion completely (see TI 
for N-body simulations in the presence of the gas drag). For this reason we believe that our results 
with no velocity evolution are applicable to the problem of planet accumulation in many cases. 

All the stirring coefficients are functions of the vertical and horizontal velocity dispersions 
of planetesimal population. Stewart & Ida (2000) have in particular shown that (Pys), (Qvs) 
s~ 2 lns for s ^ 1 [we used this result to derive equation (22)]. It means that to close properly 
the problem we need to couple equations for the velocity evolution to equation (16). In doing 
so, one should bear in mind that random motions are highly nonuniform in space, since stirring 
by the massive body is strongly localized within several Hill radii from it. Also, it is not clear 
that the velocity distribution of planetesimals can be adequately described by the Schwarzschild 
distribution, as is usually assumed for simplicity. For this reason we will not pursue this subject 
here and postpone its more detailed exploration to a future work. 

It is important to note however, that / oc s~ 2 lns and A oc s 2 /lns for s 3> 1, as follows from 
the equations (22) and (23). Thus, as A grows and the embryo heats up planetesimal population 
around it, the planetesimal "viscosity" decreases (increasing A even more through this velocity 
coupling), which facilitates gap opening. This only strengthens our conclusion that a gap must 
form when the condition A > 1 is fulfilled, even without knowing the details of the random velocity 
evolution of planetesimals. 

4.3. Effects of planetary migration. 

In §3 we studied gap opening around a massive body, assuming that the background distribu- 
tion of surface density of planetesimals is symmetric with respect to the position of the embryo. In 
this case we assumed that the embryo is fixed in radius and there is no migration at all. 

It is more than likely that in real protoplanetary disks there are significant surface density 
gradients, which could drive embryo migration. This could in principle introduce significant changes 
into our picture. Indeed, if the embryo is able to migrate quickly it may move out of the gap it 
starts to form and, thus, gap formation would be suppressed. 

Planetary migration will naturally occur in the course of embryo accumulation, since in the 
process of scattering planetesimals, the massive body exchanges its angular momentum with them, 
which leads to its migration. Indeed, planetesimals passing within Rh from the embryo get dis- 
placed by a distance ~ Rh, thus an embryo itself is displaced by ~ Rnimo/Mp). During the time 
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interval At approximately (T,o/mo)R'j I QAt planetesimals pass within embryo's Hill sphere. The 
surface number densities on both sides will likely be different by ~ (£o/mo)(-R-ff/a), although it 
is likely that as the embryo moves in some direction it plows planetesimals in front of it, leaving 
behind a depression, which would tend to oppose the migration [see Ward & Hourigan (1989) for a 
similar effect in gaseous disks]. Thus, our previous assumption about the surface density difference 
is likely to be an upper limit and the actual migration will be weaker. One can easily calculate 
that the rate of this "maximum" migration is 

1 da dh „S a 2 
— — = — ~ il — — — . 35) 

R H dt dt M c y ' 

The time it takes the embryo to migrate through a zone with a width equal to its own Hill radius 
is then r2 _1 M c /(So« 2 ) and should be longer than t open given by equation (24) if a gap is to be 
maintained. Thus, the necessary condition here is 

W c ) >mT- (36) 

If, say, Eoa 2 /M c = 1CT 5 , (and A > 1) then all bodies with masses > 10 18 g will open a gap faster 
than they migrate through it. 

Another type of migration could arise if the whole system is immersed in a massive gaseous 
disk, as should be the case in the early stages of protoplanetary evolution. In this case Goldreich 
& Tremaine (1980) demonstrated that it takes time ~ J7~ 1 (/i 2 Aa/a 3 )(M c 2 /M p S s a 2 ) for a planet 
to migrate a distance Aa in the radial direction, where h is the disk thickness, determined by its 
temperature, and T, g is the surface density of gaseous disk. In our case, the relevant lengthscale is 
again Aa ~ Rh, thus migration timescale is 

Zjgtl (1 

Comparing these two timescales we obtain: 

tmig Mp M c n 2 Vl0 24 gy 10- 3 M Q \M e J {SOJ- W 

Thus, migration due to the interaction with a gaseous disk is unlikely to have an important 
effect on gap opening. 



5. Summary 

We studied the possibility that planetary formation due to the accretion of planetesimals could 
be significantly slowed or even stopped, due to gap formation around a forming planetary embryo, 
caused by the strong gravitational perturbations of planetesimals in its vicinity. We find a critical 
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parameter A which describes the importance of gap formation [eq. (19)]. Predictions based on this 
parameter were compared with numerical N-body simulations of this process (IM93, TI97), and 
they are in good agreement. 

Only the case of a single mass distribution of the particles in a disk was studied here. But it 
is plausible that our basic results hold true even if a distribution of masses exists. Only the char- 
acteristic planetesimal mass entering this parameter should be chosen carefully, and this question 
merits further investigation [see Kokubo & Ida (1996), (1998) for some numerical results]. 

Our findings were confirmed by solving the evolution equation (16) neglecting the velocity 
dispersion evolution of planetesimals in the disk. The kinematic properties of the planetesimal 
population are important in this sort of study, and the surface density evolution and velocity 
evolution of planetesimals are closely related. Nevertheless, we hope to have grasped the main 
qualitative features of the evolution of the distribution of guiding centers even without keeping track 
of the velocity evolution. The instantaneous density of planetesimals depends on their kinematic 
properties and should experience at least a decrease by a factor of several in the vicinity of the 
embryo, leading to slowing down the accretion (and gas drag and inhomogeneous velocity evolution 
could clear out the gap completely) . We also stress that our results for a cold disk provide an upper 
limit to the embryo mass required to open a gap. 

If the disk is not uniform, migration of the embryo itself also likely to occur in the course 
of its interaction with planetesimals or with the more massive gaseous disk from which the whole 
embryo-planetesimal system originally condensed. We have estimated how this could affect our 
results and show that gap formation is likely to occur even when migration is present. 

Finally, the evolution equation itself, coupled with the equations of planetesimal velocity evo- 
lution, provide us with a powerful tool to study the formation and evolution of planetary embryos. 
Since our results seem to be in good agreement with N-body simulations, we may use this machin- 
ery for other problems of a similar nature. An obvious example is the coupled evolution of several 
planetary embryos. 
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A. Derivation of equation (6). 

Following Petit & Henon (1987b), we first calculate the current of particles with mass m\ 
(initially located at the guiding-center radius r\) through the circle of radius r, due to the inter- 
action with a single particle of mass m,2, located at guiding-center radius r 2 . The characteristic 
dimensionless relative distance for the two particles involved in this interaction 

h= J ri " r2 ! /3 - ( A1 ) 

It was demonstrated by Henon & Petit (1986) that in coordinates normalized in this way, the 
equations of relative motion of particles do not depend upon their masses. For the particle mi 
initially located to the left of the boundary at r, to cross it to the right one needs 

Ah> Ah min = ii '—j, , A2 

r m 2 (m + fi 2 ) 1/s 

where the factor (mi + m2)/m2 arises because h describes relative motion of particles. 

Then the left-to-right part of this flow of particles in the mass interval (mi, mi + dm{) during 
the time dt is obviously given by 

r oo 

J dr 1 2\A\\r 1 -r 2 \ J d(Ah)x27rnP(h,Ah)N(m 1 ,n,t)dtdm 1 . (A3) 

Here \A\\ri — r2\ is the local shear velocity between the interacting particles. 

Summing over all possible positions and masses of particles ra2 we obtain the total left-to-right 
flow of particles in the range (mi, mi + dm\): 

co oo r oo 

(AJ+) = dmi J dm 2 j dr 2 J dri j d(Ah) 

-oo -oo Ah min 

xN(m 1 ,ri,t)N(m2,r 2 ,t)P(h,Ah)4TTr 1 \A\\r 1 - r 2 \dt, (A4) 

This formula coincides with equation (28) of PH. We will further denote fi = f(rrii,...),i = 1,2 
for any function / for brevity, and replace the factor r\ under the integral with r (since r\ weakly 
varies on the scale of the Hill radius) . 

Making the change of variables from ri,r 2 to R, h given by [cf. equations (30) and (31) of PH] 

n=r + R, (A5) 
r 2 = r + R - (m + fi 2 ) 1/3 rh, (A6) 

one can reduce equation (A4) to 

oo oo oo 

(AJ+) = 4n\A\r 3 dt dm x J dm 2 (m + fi 2 ) 2/3 J dh J dR J d(Ah) 

-oo -oo Ah min 

xNi(r + R)N 2 [r + R - (m + /i 2 ) 1/3 rh}P(h, Ah)\h\. (A7) 
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We can change the order of integration over d(Ah) and dR to get for (AJ+) 



oo oo 



(AJ+) = 4ir\A\r 3 dt d mi J dm 2 (fii + fi 2 ) 2/3 J dh\h\ J d(Ah)P(h,Ah) J dR 



with 



-oo o D 

xN 1 (r + R)N 2 [r + R- (jn + fi 2 ) 1/3 rh], (A8) 



d = - < r^3 Ah - ( A9 ) 



Considering now the right-to-left flow of particles through r one can get for this component of 

flux 

oo oo D 

(AJ_) = 4ir\A\r 3 dt dm x J dm 2 (m + fi 2 ) 2/3 J dh\h\ J d(Ah)P(h,Ah) J dR 

— oo — oo 

xN 1 (r + R)N 2 [r + R - ( Pl + fi 2 ) 1/3 rh]. (A10) 



Now, the total flux of particles through the boundary at r is given by {A J) = (AJ + ) — (AJ_). 
Then the equation of evolution is obtained by setting 

^[2itrN 1 {m,r,t)]dtdm 1 = ~^^ L - (All) 

Here we do not differentiate r 3 in the right-hand side because it varies only weakly on the scale 
of Hill radius. Then the right-hand side of (All) contains d(N\N 2 )/dr which obviously equals 
d{N\N 2 ) l 'dR. Taking this into account, one can trivially perform an integration over R in the 
r.h.s. of (All) to obtain finally equation (6). In deriving it we have also taken into account that 

d(Ah)P{h,Ah) = 1. (A12) 



One should note that in deriving (A7) we assumed that (NiN^^ = (Ni)^(N 2 ) ( p, where (g)^ 
means averaging quantity g over the azimuthal angle. This might not be true in the planetary disks 
which are cold and have a large viscosity (Spahn & Sremcevic 2000), or during the initial stages 
of the gap development in planetesimal disk. However, we believe that it is unlikely to affect our 
results, since for hot disks and late times this separability assumption should be adequate. Thus, 
it is possible that our numerical results presented in §3 are somewhat different quantitatively at 
very early times from what a more detailed theory would predict. But we believe that our principal 
conclusions remain unchanged. 
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B. Form of h(ti) used in PH. 



Petit & Henon (1987b) have solved numerically the Hill equations in the case when the initial 
random motion of interacting particles is small. In this case the outcome of the interaction between 
two particles is deterministic, and they obtained the following form for the function h(h'), where h 
is the initial difference of semimajor axes of particles and h' is the final value of the same quantity, 
normalized by a[(m\ + m2)/M c ] 1//3 , where mi and rri2 are the masses of interacting bodies: 



' h + 3.34377(/i 5 + 0.2/i 4 - 3.14/i 3 )" 



h'(h) 



if h > 1.75622, 
if 1.75622 > h > 1.6777, 



350/i 2 - 1204/i + 1038.5, 
2895.903/i 5 - 20454.39/i 4 + 57671.78/i 3 

-81146.35/j 2 + 56984.57/j - 15977.97, if 1.6777 > h > 1.2219, 

-1832.5/r 2 + 4361.35/j - 2596.5, if 1.2219 > h > 1.17, 

-h - 4.107085921(/i- 1 + 4/i 4 ) exp(-5.58505361/i 3 ), if 1.17 > h > 0. 



(Bl) 



